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Résumé 

Les résultats de simulations Monte-Carlo du modèle de Hubbard à deux dimensions à demi-remplissage sont 
présentés. La méthode de simulation utilisée est celle proposée par Suzuki et al, et Hirsch et al.. Les états générés 
par cette méthode sont des états de base dans la représentation en nombres d'occupation basée sur les états de 
Wannier localisés autour de chaque site du réseau carré. Les configurations des fermions sur le réseau 2D réel peuvent 
être observées. Un facteur d'antiferromagnétisme a pu être défini et calculé pour chaque température. Les courbes 
d'énergie, de chaleur spécifique, de conductivité et de facteur d'antiferromagnétisme en fonction de la température 
sont présentées pour différentes valeurs de l'interaction répulsive coulombienne sur site U (avec t = 1). Pour les 
faibles valeurs de U le modèle présente une transition métal-isolant à basse température. Cette transition correspond 
à une transition paramagnétique-ferromagnétique du premier ordre. En effet, pour ces valeurs de l'interaction, les 
courbes d'énergie comportent un décalage qui est caractéristique d'une transition du premier ordre. De même, les 
courbes de conductivité font apparaître un phénomène d'hystérésis qui confirme la nature de la transition. Il apparaît 
un changement de comportement pour U/t ~ 3.5. La transition change de nature. Pour les valeurs U > 3.5 le 
passage ferromagnétique-paramagnétique subsiste toujours mais sans influence sur l'énergie et la chaleur spécifique. 
La transistion métal-isolant n'existe plus, la conductivité restant faible. Les réseaux d'isothermes des différentes 
grandeurs , en fonction de U /t, font apparaître une transition qui semble correspondre à la transition métal-isolant 
de Mott. Ces résultats permettent de tracer un diagramme de phase avec deux lignes de transition du premier ordre. 

Abstract 

The Quantum Monte-Carlo simulations of the two-dimensional Hubbard model are presented for the half filling. 
The method based on the direct-space proposed by Suzuki and al., and Hirsch and al. was used. The states generated 
by this method are basis states in occupation number représentation built with Wannier states localised on each site 
of the square array. The configurations of fermions can be observed on the real 2D array. An antiferromagnetic factor 
is defined and calculated for each température. The curves of energy, spécifie heat, conducivity and antiferromagnetic 
factor are presented for différent values of the répulsive coulombian on site interaction U. There is a metal-insulator 
transition at low température for small values of U. This transition corresponds with a paramagnetic-ferromagnetic 
first order transition. Indeed, for thèse interaction values, the energy curves show a gap which is a characteristic of 
a first order transition. An hystérésis phenomenon appears on the conductivity curves. There is a behaviour change 
for U/t ~ 3.5. For the values U > 3.5 there is ferromagnetic-paramagnetic change without observable effect on the 
energy and the spécifie heat. The metal-insulator transition does not exist any more, the conductivity stays very 
small. Isotherms of the physical quantities versus U /t show a transition which seems to be the metal-insulator Mott 
transition. Thèse results allow to draw a phase diagram with two first order transition fines. 

PACS numbers: 71.10.Fd Lattice fermion models, 71.27.+a Strongly correlated électron systems;heavy fermions, 
71.30.+h Metal-insulator transitions and other electronic transitions 



I Introduction 

Le modèle de Hubbard à deux dimensions est certai- 
nement le modèle le plus simple pour décrire le com- 
portement d'un système de fermions fortement corrélés. 

II semble suffisant pour expliquer la transition de phase 
métal-isolant prévue par Mott. La solution exacte de ce 
modèle n'étant pas connue il a été étudié avec de nom- 
breuses méthodes analytiques et numériques [1—17]. 

Un grand nombre de ces résultats approchés ont été obte- 
nus par des simulations numériques réalisées avec la mé- 
thode de Monte-Carlo quantique appelée méthode du dé- 
terminant. Cette méthode basée sur la formule de Trotter- 
Suzuki utilise la transformation de Hubbard-Stratonovich 
[20,21]. 

Dans cet article nous présentons les résultats de simulation 
obtenus par une autre méthode de Monte-Carlo, proposée 
par Suzuki [22,23] et Hirsch [24,25]. Cette méthode per- 
met de générer des configurations des fermions dans l'es- 
pace réel ce qui n'est pas le cas pour les autres méthodes 
numériques. Cette nouvelle méthode, qui est la version 2D 
de la méthode des u world Unes, nous a permis de retrouver 
le comportement général prévu par les autres méthodes et 
d'obtenir des résultats complémentaires. 
Cet article est organisé comme il suit : 

- dans la section 2 nous présentons les paramètres du 
modèle simulé, 

- la méthode est rappelée dans la section 3, 

- les résultats sont présentés dans la section 4, 

- nous analysons les résultats dans la section 5, 

- la conclusion est donnée dans la dernière section. 



où D est l'opérateur densité et Z est la fonction de parti- 
tion. 
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j3 = 1/ksT est l'inverse de la température. L'interac- 
tion de saut "t" introduit des éléments non-diagonaux, en 
conséquence les états propres du hamiltonien H (Eq.(JTJ), 
ne sont pas les états de base. Aussi le calcul de la trace dans 
les équations Eqs.Q et (fj| est très problématique. La mé- 
thode proposée par Suzuki et Hirsch permet de contourner 
cette difficulté [22,24]. 

L 'hamiltonien du système est décomposé en plusieurs sous- 
hamiltonicns H r . A cause des relations d'anticommuta- 
tion certains de ces sous-hamiltoniens ne commutent pas. 
Cette décomposition n'est pas totalement arbitraire, les 
sous-hamiltoniens sont choisis tels qu'ils peuvent être eux- 
mêmes décomposés en plusieurs sous-systèmes sans site 
commun dont les hamiltoniens K r ^ commutent. 

p 

H = J2 H r [H r ,H r ,]^0 r^r 1 (5) 

r=l 

H r = J2 Kr > k [K r , k ,K rtk ,]=0 Vk,k' (6) 
fc=i 

Les décompositions précédentes sont réalisées dans l'in- 
tention d'utiliser la formule de Trotter, Eq.f7|), qui permet 
de contourner le problème de la non-commutativité des 
sous-hamiltoniens H r . 



2 Modèle de Hubbard 

Le modèle étudié est un réseau carré. Les états mono- 
particulaires utilisés pour construire les états de base 
dans la représentation en nombre d'occupation sont les 
états de Wannier localisés sur les sites. L 'hamiltonien de 
Hubbard est 
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En utilisant la formule de Trotter, la fonction de partition 
Z s'écrit 

Z = lim Z n (8) 



où Z n est un approximant de la fonction de partition 
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où les opérateurs c\ a et a, a sont les opérateurs de création En insérant np ensembles complets d'états de base entre 
et d'annihilation d'un fermion de spin a sur le site i. (i, j) les opérateurs, l'approximant Z n devient 



indique que la sommation porte sur les premiers voisins. 



3 Méthode numérique 

La méthode de simulation est celle développée dans la 
référence [26]. Son principe est présenté dans les références 
[23-25]. Nous en rappelons les principales caractéristiques. 

3.1 Principe de la méthode de simulation 

Dans l'ensemble canonique la valeur moyenne d'une ob- 
servable, O, est donnée par 



Zn = CXP (~~ H P ) \*np-l) 
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[\& Q ] représente la configuration des np états et peut 
être considéré comme l'état d'un système classique de di- 
mension (d + 1), où d est la dimension du système quan- 
tique étudié, {[^q]} indique que la somme porte sur toutes 
les configurations possibles. Ceci est équivalent à diviser 
le temps imaginaire t, tel que < t < 0, en n intervalles 
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de durée Ar = (3/n. On détermine un approximant de 
l'énergie du système à partir de la relation 



En remplaçant Z n par Eq. (|10[) . on obtient : 
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L'indice r du sous-hamiltonien qui apparaît dans chaque 
élément de matrice est fonction de l'indice de l'état j, tel 
que r = 1 + ( j mod p) . Le calcul de la trace impose des 
conditions périodiques sur les états, ainsi le ket |^o) cor- 
respond au ket l^np). Les facteurs P n vérifient 
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On peut ainsi considérer que chaque configuration [\I> a ] du 
système de dimension (d+ 1) a une énergie £7 n et 
un facteur de probabilité P n ([^ a ]). Le calcul de la valeur 
moyenne de l'énergie du système est réalisé en utilisant une 
méthode de Monté-Carlo telle que l'algorithme de Metro- 
polis. 

En utilisant la relation ©, la valeur de U n s'écrit : 



où l'état \tfj,i) est l'état du sous-système d'indice l du sous- 
hamiltonien d'indice r de l'intervalle de temps d'indice i. 
Le symbole — > signifie que l'expression du membre de 
gauche est remplacée par celle du membre de droite. Cette 
opération implique que l'espace des états du système est 
considéré comme le produit tensoriel des espaces des états 
de tous les sous-sytèmes. Ceci n'est pas correcte pour un 
système de fermions, en effet, à cause de l'antisymétrie des 
états, les opérateurs création et annihilation sont définis 
dans l'espace des états de tout le système. Il n'y a pas une 
égalité stricte entre ces deux expressions. Implicitement 
ceci signifie que les nombres d'occupation des autres états 
monoparticulaires ne sont pas pris en compte. 
Etant donné la forme particulière des hamiltoniens, cette 
approximation peut, selon les valeurs des nombres d'oc- 
cupation des états monoparticulaires, changer le signe de 
certains éléments non-diagonaux des K r j. 
Si le nombre de types de sous-systèmes est faible, cette 
simplification permet de réduire considérablement les cal- 
culs en diminuant les dimensions de l'espace des états où 
sont mené ces calculs. En effet, la diagonalisation du ha- 
miltonien H de tout le système dans l'espace des états 
de dimension 2 N , où N est le nombre de sites de tout le 
système, est remplacée par la diagonalisation des hamilto- 
niens K ry k de quelques sous-systèmes, dont les espaces des 
états sont de très petites dimensions comparées à 2 N . 
On applique la même factorisation pour le numérateur de 
l'équation Eq. (fTT|) . Les expressions calculées sont finale- 
ment : 
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où N p est le nombre de configurations retenues, a est l'in- 
dice des configurations qui remplace [ty a ] et i est l'indice 
des intervalles de temps. Les indices i, j et r sont reliés par 
le principe de numérotation des np états et vérifient 
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La décomposition des sous-hamiltoniens H r en sous- 
systèmes permet de réaliser une simplification importante. 
Cette simplification nécessite une approximation. Tous 
les hamiltoniens K r j des sous-systèmes d'un même sous- 
hamiltonien H r commutent et chaque hamiltonien K r> i 
agit seulement sur l'état d'un seul sous-système aussi 
chaque dénominateur de l'expression (fTT|) peut s'écrire 



Certains facteurs de probabilité P n ([*«]) où P' na sont 
négatifs : c'est le "problème du signe". Nous avons montré 
dans la référence [26] qu'on peut ne pas tenir compte de 
ce signe. 

La chaleur spécifique est déterminée à partir de la formule 



(24) 
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où ks est la constante de Boltzman. 

En utilisant une méthode similaire à celle utilisée pour U' n 
on obtient une formule qui permet de calculer la chaleur 
spécifique à partir des fluctuations 
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où C est la chaleur spécifique molaire et R la constante 
des gaz parfaits. 

3.2 Conductivité statique 

La méthode de simulation ne permet pas de faire une 
détermination rigoureuse de la conductivité statique. En 
effet, l'opérateur densité de courant défini par Scalapino 
and al. [28] n'est pas diagonal et ne peut pas être dé- 
composé en opérateur relatif à chaque sous-système [25]. 
Cependant, on obtient une valeur approchée de la densité 
de courant au temps imaginaire r = ^j, en calculant : 

jx = ^ n^i-j (1 - n a> ir.j) n a ,i>;j> 

a 

(1 - n cr,l;j') - n a ,i'.j (1 - n a ,i ; j) 



(26) 



j' = j +p où p est le nombre de sous-hamiltoniens. Les co- 
ordonnées des sites l et V vérifient xy = xi + 1 et yy = yi. 
rio-j-j est le nombre d'occupation de l'état monoparticu- 
laire \a, l) du réseau 2D d'indice j. Cette formule de calcul 
de j x correspond à l'expression 



(27) 



où l'état \tpj) correspond à l'état des sites letl'. La conduc- 
tivité calculée est : 

j np JV 

A** (m) = — yy£ j 3x {l-,3)3x{l-,3+m) (28) 
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où N est le nombre de sites du système. La conductivité 
statique est obtenue après une transformation de Fourier 
discrète. 

3.3 Facteur d'antiferromagnétisme 

Pour chaque température la fonction de corrélation de 
spin c(l Xl l y ) est calculée 



c {lx,l v ) = ((K,t - «a) (m+i, r - n i+Li))) 



(29) 



(• • • ) signifie que, pour chaque configuration retenue, on 
calcule la moyenne pour tous les N sites de chaque ré- 
seau 2D et pour tous les réseaux 2D de cette configura- 
tion. La fonction finale correspond à la moyenne sur toutes 
les configurations retenues à une même température. I re- 
présente le changement d'indice des sites provoqué par la 
translation de vecteur (l x , l y ) sur le réseau. Le facteur d'an- 
tiferromagnétisme (paramètre d'ordre de Néel) est déduit 
de la fonction de corrélation de spin tel que 



(Ix+ly) 



(30) 



Dans le cas d'un ordre antiferromagnétique ce facteur est 
tel que Ja = 1, par contre, il n'est pas rigoureusement 
nul pour une absence totale d'ordre. En effet, le terme 
correspondant à l x = l y = n'est pas éliminé dans notre 
définition de ce facteur. Cette contribution est de l'ordre 
de - l/N. 
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FiG. 1 - (Couleur en ligne) Décomposition du réseau 
carré en deux sous-hamiltoniens et en sous-systèmes avec 
les conditions aux limites périodiques. Les sous-systèmes 
gris (rouge) appartiennent au sous-hamiltonien 1, les sous- 
systèmes noirs appartiennent au sous-hamiltonien 2. 

4 Résultats 

4.1 Paramètres des simulations 

Le modèle étudié est un réseau carré contenant N — 
6x6 sites avec les conditions aux limites périodiques. 
Tous les sous-systèmes sont identiques et sont composés 
de quatre sites. Ainsi la méthode ne nécessite la diagona- 
lisation que d'une seule matrice de dimension 16 x 16. Ces 
sous-systèmes sont regroupés dans deux sous-hamiltoniens 
(p = 2). La figure [T] montre la décomposition du système 
en sous-systèmes. Toutes les simulations ont été réalisées 
avec huit intervalles de temps imaginaire (n — 8). L'in- 
teraction de saut est fixée à t = 1. Les premiers voisins 
d'un site sont les quatre sites correspondant aux trans- 
lations (Ai = ±1, Ay = 0) et (Ax = 0, Ay = ±1). Les si- 
mulations ont été programmées pour des valeurs de l'inter- 
action répulsive coulombienne variant de U = à U = 8. 
La majorité des simulations ont été réalisées pour le demi- 
remplissage, c'est à dire dix-huit spins up et dix-huit spins 
down (18 t +18 |). Des simulations, pour U = 0, ont été 
menées pour un seul spin up (1 t +0 J.) et dix- huit spins 
up sans spin down (18 | +0 |) pour comparer les résultats 
afin de vérifier la cohérence de la méthode. La plupart 
des simulations ont consisté en neuf cycles de descente- 
montée de température, avec 100 points de température 
en progression géométrique pour chaque montée ou des- 
cente de température. Pour chaque grandeur étudiée trois 
moyennes ont été calculées : la moyenne en montée, la 
moyenne en descente et la moyenne totale. Ceci a permis 
de mettre en évidence un faible phénomène d'hystérésis 
entre la montée et la descente de température dans une 
certaine plage de température. 

4.2 Test de la méthode (U=0) 

La figure [2] montre les courbes d'énergie et de chaleur 
spécifique par spin en fonction de la température pour 
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Température 

FiG. 2 - (Couleur en ligne) Courbes d'énergie et de chaleur 
spécifique par spin en fonction de la température pour U = 
0. 




Température (kT) 

FiG. 3 - (Couleur en ligne) Courbes d'énergie par site en 
fonction de la température pour différentes valeurs de U. 



FiG. 4 - (Couleur en ligne) Courbes de chaleur spécifique 
molaire en fonction de la température pour différentes va- 
leurs de U. 

U = et les trois remplissages suivants : (1 | +0 j), 
(18 î +0 4) et (18 î +18 |). Pour U = le demi- 
remplissage (18 t +18 |) correspond à deux gaz quan- 
tiques non-corrélés qui ont les mêmes énergies et chaleurs 
spécifiques. Ceci explique pourquoi les courbes sont iden- 
tiques pour les remplissages (18 t +0 J.) et (18 î +18 |). 
L'énergie et la chaleur spécifique par un seul spin sur le 
réseau (1 t +0 |) sont les doubles des valeurs obtenues 
pour les remplissages précédents. Le niveau du fondamen- 
tal est égal à — 4t comme il est prévu. En effet cet état 
correspond à : 

1 N 

l*o) = -7= l ni 'T = °' ■ ■ ■ ' = n ^ = °> 

v i—l 

«i+i,î = 0, . . .,njv,i = 0) (31) 
dont l'énergie est donnée par 

ff|* > = -t£c} iîCi>T |¥o> = -4t|¥o) (32) 

Le rapport deux entre les grandeurs calculées pour le 
remplissage (1 | +0 J.) et les deux autres remplissages 
((18 î +0 |) et (18 î +18 I)) peut être justifié en considé- 
rant que pour ces deux derniers remplissages chaque spin 
ne peut se déplacer, en moyenne, que sur deux sites voisins 
alors que dans le cas d'un seul spin sur le réseau le spin a 
quatre possibilités de saut. 

Il est intéressant de remarquer que dans le cas d'un seul 
spin sur le réseau le problème du signe n'existe pas, de 
même, l'approximation relative aux nombres d'occupation 
non pris en compte dans la relation (|19p n'est pas utilisée. 
Dans ce cas la méthode n'utilise aucune approximation ex- 
cepté le découpage du temps imaginaire lié à la formule de 
Trotter. La bonne concordance de ces résultats confirme 
la cohérence de la méthode de simulation. 

4.3 Résultats des simulations 

Les figures El Hl [5l [6] et [7] montrent les courbes moyennes 
d'énergie, de chaleur spécifique, de conductivité, du facteur 
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FlG. 5 - (Couleur en ligne) Courbes de conductivité en 
fonction de la température pour différentes valeurs de U. 
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FlG. 6 - (Couleur en ligne) Courbes du facteur d'antiferro- 
magnétisme en fonction de la température pour différentes 
valeurs de U. 




0,001 



0,01 



0,1 i 

Température (kT) 



100 



FlG. 7 - (Couleur en ligne) Courbes de double ocupation 
en fonction de la température pour différentes valeurs de 
U. 



d'antiferromagnétisme et de double occupation (fraction 

de sites occupés par une paire de spins (ÎJ.)). Ces courbes 

correspondent aux moyennes totales sur les montées et les 

descentes des neuf cycles de température. 

Pour U = 0, le facteur d'antiferromagnétisme est j a ~ 

1/JV, c'est à dire qu'il n'y a pas d'ordre magnétique quelle 

que soit la température. 

Le facteur de double occupation 



D = J^( n *A n i,l) 



(33) 



peut être relié à la valeur moyenne du carré de l'aiman- 
tation locale (rn^). En effet, si le moment magnétique de 
chaque particule est choisi tel que s z = ±1 alors, pour 
chaque réseau 2D : 
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v 
(34) 

où 7V| = ^2 ?ii,î et Ni = ^2 n%,i sont les nombres de spins 
up et down, n p = ^2 rijj n,- j_ est le nombre de paires dans 
un réseau 2D. Dans le cas présent, où iVj = TYj, = N/2 = 
18, on obtient 

(ml) =1-2D (35) 

On remarque qu'à basse température (kT < t), pour U = 
0, le nombre de paires est exactement n p = 9, c'est à dire 
que les spins de chaque orientation sont distribués de façon 
totalement aléatoire sur le réseau 2D. Ceci est en parfait 
accord avec le résultat relatif au facteur d'antiferromagné- 
tisme. 

On remarque un point de croisement, pour kT ~ 2.7t, sur 
le réseau des courbes de chaleur spécifique de la figure [U 
Les figures 151 à ITT1 représentent les courbes de chaleur spé- 
cifique, de conductivité, de facteur d'antiferromagnétisme 
et de carré de l'aimantation locale pour les basses tempé- 
ratures avec une échelle de température linéaire. Pour les 
faibles valeurs de U, les parties décroissantes des courbes 
de conductivité sont assez bien approximées par des fonc- 
tions de la forme 



a(kT) 



a-b (kT) 
(kTf+c 



(36) 



où d est petit. 



4.4 Discontinuité des courbes d'énergie 

Pour les températures inférieures à environ kT ft = 3.5 
les courbes d'énergie présentent un léger décrochement qui 
correspond au pic de chaleur spécifique. Des simulations 
avec 200 points de température dans des plages de tem- 
pératures autour ce ces décrochements ont été réalisées 
pour U = 0M, U = lt, U = 2t, U = M, U = 3.3t, 
U = 3.5t et U = 4t. Les courbes d'énergie correspon- 
dantes (qui ne sont pas des courbes moyennes) sont pré- 
sentées dans les figures [Ï2] à [TTl Les courbes moyennes de 
la figure [3] sont reproduites pour comparaison. Les barres 
d'erreur de ces courbes permettent d'évaluer les fluctua- 
tions. On remarque, sur toutes les courbes, que l'ampli- 
tude de ces barres d'erreur est faible en dehors de petites 
plages de température. Ces plages de température corres- 
pondent à des zones de fluctuation sur les courbes avec 200 
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FiG. 8 - (Couleur en ligne) Courbes de la chaleur spéci- 
fique pour les basses températures. 
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FiG. 9 - (Couleur en ligne) Courbes de la conductivité 
pour les basses températures. 
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FiG. 10 - (Couleur en ligne) Courbes du facteur d'antifer- 
romagnétismc pour les basses températures. 
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FiG. 12 - (Couleur en ligne) Courbes d'énergie pour U = 
0.5t. La courbe en tiretés noirs correspond à la courbe 
moyenne de la figure [3] avec ses barres d'erreur. 



points de mesure. On constate que dans ces zones, les éner- 
gies fluctuent entre deux valeurs limites qui correspondent 
aux extrapolations des portions de courbe sans fluctuation. 
L'amplitude des fluctuations, qui correspond approximati- 
vement, au décrochement augmente de U = 0.5t kU = 3t, 
puis diminue rapidement quand U augmente. Il n'y a plus 
de fluctuation pour U — 3.5t. L'amplitude maximum, es- 
timée pour U = 3i, est d'environ 44 10~ 3 t. Pour U — 3.3i 
la valeur des fluctuations d'énergie est d'environ 14 10~ 3 i, 
et l'extrapolation des courbes est trop "hasardeuse" pour 
présenter un intérêt. 

4.5 Phénomène d'hystérésis 

Les figures [TÏÏ] à [53] présentent les courbes de conducti- 
vité moyennes calculées en descente et en montée de tem- 
pérature. Ces courbes moyennes font apparaître un phé- 
nomène d'hystérésis entre la montée et la descente de tem- 
pérature dans la zone de transition, pour les valeurs de U 
inférieures à 3.5t. Bien qu'il s'agisse de courbes moyennes 
et que les courbes en montée de température soient com- 
prises dans les barres d'erreur des courbes tracées en des- 
cente, la répétition systématique de ce phénomène pour 
toutes les valeurs de U < 3.5t indique que ce n'est pas 
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FiG. 13 - (Couleur en ligne) Courbes d'énergie poui 
U = lt. La courbe en tiretés noirs correspond à la courbf 
moyenne de la figure OU avec ses barres d'erreur. 
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FiG. 14 - (Couleur en ligne) Courbes d'énergie poui 
U = 2t. La courbe en tiretés noirs correspond à la courb< 
moyenne de la figure OU avec ses barres d'erreur. 
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FiG. 15 - (Couleur en ligne) Courbes d'énergie p< 
U = 3t. La courbe en tiretés noirs correspond à la coui 
moyenne de la figure [3] avec ses barres d'erreur. 
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FiG. 16 - (Couleur en ligne) Courbes d'énergie pour 
U = 3.3t. Les barres d'erreur correspondent à la courbe 
moyenne de la figure 03 
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FiG. 17 - (Couleur en ligne) Courbes d'énergie pour 
U = 3.5t. Les barres d'erreur correspondent à la courbe 
moyenne de la figure [3J 
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FiG. 18 - (Couleur en ligne) Courbes de conductivité en 
descente et en montée de température pour U — 0.5t. Les 
barres d'erreur correspondent à la courbe en descente de 
température. 
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Température (kT) 

FiG. 19 - (Couleur en ligne) Courbes de conductivité 
descente et en montée de température pour U = Xt. 1 
barres d'erreur correspondent à la courbe en descente 
température. 




Température (kT) 

FiG. 20 - (Couleur en ligne) Courbes de conductivité 
descente et en montée de température pour U = 2t. 1 
barres d'erreur correspondent à la courbe en descente 
température. 



un phénomène aléatoire. Cette mise en évidence d'un p 
nomène d'hystérésis doit, cependant, être considérée a 
prudence. En effet, conformément au principe des r 
thodes de Monte-Carlo, les états générés à une tempé 
ture donnée ne doivent pas être fonction des états pré 
dents. Ce phénomène d'hystérésis apparaît également 
les courbes du facteur d'antiferromagnétisme. 

4.6 Résultats en fonction de U/t 

Comme il a été indiqué précédemment, les simulatic 
ont été réalisées en programmant des cycles de tempéri 
ture. Les résultats obtenus ont permis de tracer les v; 
riations de différentes grandeurs en fonction de la tempi 
rature. Ceci est parfaitement adapté pour mettre en év 
dence des transitions de phases induites par les variatioi 
de température. Par contre, la transition de Mott attei 
due doit être, principalement, produite par les variatioi 
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FiG. 21 - (Couleur en ligne) Courbes de conductivité en 
descente et en montée de température pour U — 3t. Les 
barres d'erreur correspondent à la courbe en descente de 
température. 
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FiG. 22 - (Couleur en ligne) Courbes de conductivité en 
descente et en montée de température pour U = 3.3t. Les 
barres d'erreur correspondent à la courbe en descente de 
température. 




Température (kT) 

FiG. 23 - (Couleur en ligne) Courbes de conductivité en 
descente et en montée de température pour U = 3.5t. Les 
barres d'erreur correspondent à la courbe en descente de 
température. 
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U/t 

FiG. 24 - (Couleur en ligne) Isothermes d'énergie en fon< 
tion de U/t. Seuls les points sur l'isotherme kT = 0.210 
ont été reportés. Les points sur les autres courbes corres 
pondent aux mêmes valeurs de U/t. 

de l'interaction coulombienne U/t. Il est donc souhaitât»' 
de tracer l'évolution des différentes grandeurs en fonctio 
de U/t. Ceci est indispensable si la ligne de transition et 
parallèle à l'axe des températures dans le plan (U/t, kT/t 
c'est à dire si elle n'est pas traversée en programmant de 
cycles de températures. Les résultats des simulations or 
donc été traités pour obtenir les courbes souhaitées. Le 
isothermes d'énergie , de chaleur spécifique, de conduct 
vité, du facteur d'antiferromagnétisme et de double occi 
pation sont présentées dans les figures l2"4l à I5U1 




FiG. 25 - (Couleur en ligne) Isothermes de chaleur spé- 
cifique en fonction de U/t. Pour 0.2656 > kT > 0.13t les 
courbes se superposent dans la zone de transition. 



5 Discussion 

5.1 Transition métal-isolant du premie 
ordre 

Les courbes de chaleur spécifique de la figure [U ci 
conductivité des figures [S] et [5] et du facteur d'antiferromf 
gnétisme des figures [ïïj et [TUJ font clairement apparaître un 
transition de phase, quand la température varie, pour le 
valeurs de l'interaction coulombienne telles que U < 3.5 
La phase basse température est une phase isolante ar 
tiferromagnétique. La phase haute température est un 
phase conductrice sans ordre magnétique. Les courbes d 
conductivité à basse température de la figure [5] montrer 
que cette phase a un comportement métallique. La trai 
sition est donc une transition métal-isolant concomittant 
à une transition paramagnétique-antiferromagnétique. Le 
courbes moyennes des figures [3] à [7] ne permettent pas d 
déterminer avec certitude la nature de cette transitioi 
Par contre, les figures fT2l à [TTl montrent clairement que le 
courbes d'énergie moyenne comportent un décrochcmcr 
caractéristique des transitions du premier ordre. Ce de 
crochement est confirmé par les fluctuations des courbe 
avec 200 points de mesure des figures [J3] à [TU En effet, ce 
courbes s'interprètent en considérant que, dans les plage 
des températures de transition, les états du système coi 
respondent à deux puits de potentiels. Ainsi, les états su< 




FiG. 26 - (Couleur en ligne) Isothermes de chaleur spé- 
cifique en fonction de U /t pour les hautes températures 
(kT > 0.2656t). 
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FiG. 27 - (Couleur en ligne) Isothermes de conductr 
en fonction de U/t. Seuls les points sur l'isotherme kl 
0.2105 ont été reportés. Les points sur les autres cour 
correspondent aux mêmes valeurs de U/t. 




Facteur d' antiferromagnétisme 



FiG. 29 - (Couleur en ligne) Isothermes du facteur d'an- 
tiferromagnétisme en fonction de U/t. 




FiG. 28 - (Couleur en ligne) Isothermes de conducth 
en fonction de U/t pour les hautes températures. 




Double occupation 



FiG. 30 - (Couleur en ligne) Isothermes de double occu- 
pation D en fonction de U/t. 
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cessifs générés par l'algorithme peuvent rester "bloqués" 
longtemps dans un des puits de potentiel. Ceci justifie que 
les valeurs de l'énergie fluctuent entre deux valeurs limites. 
Il y a donc un phénomène de métastabilité entre deux états 
macroscopiques possibles du système. Ce phénomène de 
métastabilité et le phénomène d'hystérésis observé sur les 
courbes de conductivité sont également des propriétés ca- 
ractéristiques des transitions du premier ordre. On peut 
donc conclure, sans aucun doute, que la transition est du 
premier ordre. On en déduit que les pics de chaleur spéci- 
fique ne sont pas correctes. De même, les valeurs moyennes 
des autres grandeurs ne sont pas exactes dans les plages 
de températures où il y a le phénomène de métastabilité. 
En effet, dans la zone de transition, les configurations gé- 
nérées par l'algorithme qui sont retenues pour calculer les 
différentes grandeurs, sont des états microscopiques des 
deux phases. C'est à dire que le système "oscille" entre les 
deux phases. Ce phénomène de métastabilité signifie que le 
système n'est pas à l'équilibre thermodynamique, donc la 
relation [25] qui permet de calculer la chaleur spécifique ne 
peut pas être utilisée. La singularité mathématique atten- 
due sur la courbe de chaleur spécifique d'une transition du 
premier ordre est un saut ou décrochement comme pour 
l'énergie. 

5.2 Point critique 

Comme cela a été remarqué précédemment, un change- 
ment de comportement se produit pour U ~ 3.5t. Cette 
valeur est la valeur limite d'une ligne de transition du pre- 
mier ordre dans le plan (U/t,kT/t). La température cor- 
respondante, à mi- hauteur du pic, est environ kT/t ~ 0.13. 
Ces deux valeurs sont les coordonnées du point critique 
qui termine la ligne de transition du premier ordre. Les 
températures des points de cette ligne de transition sont 
également déterminées à mi-hauteur des pics. Le réseau 
d'isothermes de la figure [2"ÏÏ1 confirme qu'il s'agit bien d'un 
point critique analogue à celui observé pour les transitons 
liquide-gaz. 

On remarque, sur les figures [53] et [3D] que les paliers des 
isothermes telles que kT/t < 0.1322 ne sont pas parallèles 
à l'axe des "x". Ceci est dû à la méthode de Monte-Carlo. 
En effet, comme cela a été indiqué précédemment, dans 
une zone de transition du premier ordre cette méthode gé- 
nère des états des deux phases possibles, ainsi les moyennes 
obtenues ne correspondent pas à une seule phase. On ob- 
tient alors des valeurs intermédiaires entre les valeurs rela- 
tives à chaque phase. Il faudrait, pour contourner ce pro- 
blème, sélectionner les configurations retenues en ne gar- 
dant que celles correspondant à une seule phase, c'est à 
dire à un puits de potentiel. Malheureusement, il apparaît 
difficile de déterminer à quelle phase correspond une confi- 
guration générée. Aussi, quels que soient les nombres de 
points reportés sur les courbes, les paliers ne seront jamais 
parfaitement parallèles à l'axe des "x". 

5.3 Transition de Mott 

Alors que les réseaux d'isothermes du facteur d'antifer- 
romagnétisme de la figure [221 et du facteur de double occu- 



pation de la figure [3U] confirment l'existence du point cri- 
tique (U/t ~ 3.5, kT/t ~ 0.13), les isothermes de conduc- 
tivité et de chaleur spécifique des figures [27] et [25] font ap- 
paître un saut de conductivité associé à un saut de chaleur 
spécifique pour les températures 0.2656i > kT > 0.13t. Ce 
saut se produit pour la valeur constante U/t ~ 3.5, à mi- 
hauteur du décrochement. Les figures l2"6l et |2"51 montrent 
que pour les températures kT > 0.2656 ce saut disparaît, 
les courbes présentent alors un point d'inflexion. Cette dis- 
continuité doit correspondre à une transition métal-isolant 
de Mott. La présence de cette transition est également ob- 
servée sur les isothermes d'énergie de la figure[5U En effet, 
ces courbes présentent un point anguleux pour la valeur de 
l'interaction coulombicnnc U/t ~ 3.5. Cependant, les ré- 
sultats ne permettent pas de déterminer de façon certaine 
la nature de cette transition. Etant donné le changement 
de comportement observé sur les courbes des figures [2"ol et 
[51] pour les températures kT > 0.2656i, on peut supposer 
que la transition est encore du premier ordre dans l'inter- 
valle de température 0.1322 — 0.2656. Le point (U/t = 3.5, 
kT/t — 0.2656) serait alors un point critique sur la ligne de 
transition de Mott. Dans cette hypothèse, les deux lignes 
de transition du premier ordre doivent se rejoindre en un 
point triple proche du point critique de la transition iso- 
lant antiferromagnétique-conducteur paramagnétique. 

5.4 Diagramme de phase 

On peut considérer que la valeur moitié de la conducti- 
vité maximum est la valeur limite entre le comportement 
conducteur et le comportement mauvais conducteur. Ainsi 
pour des valeurs de la conductivité légèrement supérieures 
le système se comporte comme un conducteur tandis que 
pour des valeurs légèrement inférieures le système peut 
être qualifié de mauvais conducteur. La ligne en pointillés 
(rouge) de la figure l3"Tl correspond à cette limite de com- 
portement. 

Les courbes de conductivité montrent que le système est 
isolant et paramagnétique, quel que soit U, quand kT > 
2t. Cependant, les courbes de la figure [7] montrent que ces 
états paramagnétiques sont caractérisés par des nombres 
de double occupation différents. De même on constate 
que le système est isolant, quelle que soit la tempéra- 
ture pour U > 8£, mais le système est antiferromagnétique 
à basse température et devient progressivement parama- 
gnétique quand la température augmente. Ce passage iso- 
lant antiferromagnétique-isolant paramagnétique se pro- 
duit sans qu'il y ait de transition, au sens thermodyna- 
mique, détectable sur toutes les grandeurs calculées dans 
nos simulations. Ceci est peut être dû à la petite taille du 
réseau. 

Comme pour le changement de comportement conducteur- 
mauvais conducteur, on peut considérer que le système 
présente une tendance au ferromagnétisme pour Ja > 0.5 
et au paramagnétisme pour /a < 0.5. Ceci permet de tra- 
cer une ligne de "cross-over" entre les deux phases magné- 
tiques. 

Selon les hypothèses précédentes le diagramme de phase 
dans le plan (U/t, kT/t) comprend, donc, une première 
ligne de transition du premier ordre pour kT/t < 0.13 et 
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FlG. 31 - (Couleur en ligne) Allure générale du diagramme 
de phase du modèle de Hubbard 2D à demi-remplissage. 
Les lignes en pointillés (rouge) et (vert) sont des lignes de 
"cross-over". La ligne en trait plein (noir) est la ligne de 
transition du premier ordre isolant antiferromagnétique- 
conducteur paramagnétique. La ligne en tiretés (bleu) cor- 
respond à la transition de Mott. 



U/t < 3.5 et une deuxième ligne de transition du premier 
ordre pour 0.2656 > kT > 0.13, cette ligne étant parallèle 
à l'axe des températures. Ces deux lignes se rejoignent en 
un point triple. Les deux lignes de "cross-over" définies 
précédemment peuvent compléter le diagramme. Toutes 
ces considérations permettent de tracer le diagramme de 
phase de la figure I3T1 

5.5 Comparaison avec quelques résultats 
de la littérature 

La transition basse température métal paramagnétique 
- isolant antiferromagnétique a effectivement été observée 
expérimentalement sur un certain nombre de composés, 
notamment V2O3 [29,31,32]. Cette transition est effective- 
ment du premier ordre. Les résultats obtenus en faisant va- 
rier la pression mettent en évidence une transition métal- 
isolant à haute température qui correspond à la transition 
de Mott. Cette transition est également une transition du 
premier ordre. Les résultats de nos simulations sont donc 
en bon accord qualitatif avec l'expérience. 
Les diagrammes de phases expérimentaux font également 
apparaître une ligne de transition entre la phase isolante 
antiferromagnétique et la phase isolante paramagnétique. 
Cette transition n'a pas été détectée dans nos simulations. 
Ces différences peuvent être, éventuellement, justifiées en 
considérant que le modèle étudié est un modèle 2D alors 
que les résultats expérimentaux concernent des matériaux 
3D. De même, l'influence de l'interaction coulombienne U 
n'est certainement pas totalement équivalente à celle de la 
pression. 

De nombreux résultats de simulations et d'études théo- 
riques du modèle de Hubbard 2D ont été publiés. Nous 
ne considérerons que quelques un de ces articles qui pré- 
sentent les courbes des mêmes grandeurs que celles calcu- 
lées dans nos simulations. De façon générale, on remarque 
que l'allure de notre diagramme de phases est proche de 
celle des diagrammes obtenus par d'autres méthodes, pour 



des modèles non frustrés. En particulier, on constate que 
le point critique que nous avons déterminé correspond au 
point triple du diagramme de phases de la référence [3]. 
Les courbes d'énergie et de chaleurs spécifiques ont égale- 
ment une allure générale très proches de celles présentées 
dans les références [4,6,8,11,14]. Mais les résultats de 
ces références ne font pas apparaître les sauts de l'énergie 
qui permettent de déterminer la nature de la transition. 
On constate également que les courbes du carré de l'ai- 
mantation locale (m^) de la figure [TT] ont la même allure 
que celles présentées dans la référence [14]. Cependant, 
une différence apparaît. En effet nos courbes présentent, à 
basse température, des variations qui correspondent à la 
transition métal-isolant du premier ordre alors que cette 
transition n'est pas détectée dans la référence. 
Grôber et al. [18] ont réalisé des simulations de Monté- 
Carlo quantique du modèle 2D de Hubbard (8 x 8) en 
faisant varier la valeur de l'interaction coulombienne U 
pour la température T = 0.33t. Ils ont mis en évidence 
une transition métal paramagnétique-isolant pour la va- 
leur critique U c sa 4i. Cette valeur critique est proche de 
la valeur U ~ 3.5i que nous avons obtenue. 



6 Conclusion 

Comme toute méthode numérique, la méthode employée 
pour nos simulations utilise des approximations qui per- 
mettent de réduire considérablement les calculs. La vali- 
dité des résultats obtenus peut être mise en doute à cause 
de ces approximations. Mais la cohérence des résultats et 
le bon accord qualitatif avec un certain nombre de résul- 
tats obtenus avec d'autres méthodes numériques, avec les 
prédictions théoriques et les résultats expérimentaux per- 
met de considérer que ces résultats sont significatifs. 
Le Modèle de Hubbard 2D étudié fait apparaître deux 
transitions : une transition du premier ordre métal para- 
magnétique - isolant antiferromagnétique à basse tempéra- 
ture et la transition de Mott métal paramagnétique-isolant 
paramagnétique pour 0.2656 > kT /t > 0.13. Par contre, 
le passage isolant antiferromagnétique-isolant paramagné- 
tique se produit sans transition de phase détectable sur 
les grandeurs étudiées. Il est possible que l'effet de taille 
ne permette pas d'observer tous les "détails" du compor- 
tement du modèle qui semble être compliqué autour du 
point critique. 
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